Formation of "Blanets" from Dust Grains around the Supermassive Black Holes in Galaxies

There is enough evidence suggesting that planets are formed in the circumstellar disks around stars. However, stars might not be the only site for planet formation. Recently, in Wada et al. (2019, hereafter Paper I), we claimed a new class of "planets" ⁵ that orbit around supermassive black holes (SMBHs) in galactic centers. Paper I theoretically investigated the growth processes of planets, from sub-micron-sized icy dust monomers to Earth-sized bodies outside the snowline in a circumnuclear disk around an SMBH, typically located several parsecs from the SMBHs. As is the case in a protostellar disk, in the early phase of the dust evolution, low-velocity collisions between dust particles promote sticking; therefore, the internal density of the dust aggregates decreases with growth (Okuzumi et al. 2012; Kataoka et al. 2013). When the size of porous dust aggregates reaches 0.1–1 cm, the collisional and the gas-drag compression become effective, and as a result, the internal density stops decreasing. Once 10–100 m sized aggregates are formed, they decouple from gas turbulence, and as a result, the aggregate layer becomes gravitationally unstable (Michikoshi & Kokubo 2016, 2017), leading to the formation of "planets" due to the fragmentation of the layer, with ten times the mass of the earth. The objects orbit the SMBHs with an orbital time of 10⁵–10⁶ yr. To distinguish them from standard planets, we hereafter call these hypothetical astronomical objects blanets. ⁶

The results reported in Paper I, however, have two major limitations. One is that the collisional velocity between the dust aggregates might become too large (>several 100 m s⁻¹ at the Stokes parameter, S_t ∼ 1). And if the collisional velocity is that large, rather than growing, the aggregates might get destroyed. In Paper I, we used the numerical experiments conducted by Wada et al. (2009). ⁷ on the collisions between the dust aggregates, wherein the critical collisional velocity (v_crit) scales with the mass m_d of the dust aggregates, as ${v}_{\mathrm{crit}}\propto {m}_{d}^{1/4}$. However, this is correct only for the head-on collisions, as stated in the paper. Moreover, Wada et al. (2009, 2013) showed that the growth efficiency of the dust aggregates depends on the impact parameter of the collisions, and as a result, v_crit does not strongly depend on the mass of the dust aggregates, if offset collisions are taken into account. They concluded that v_crit ≃ 80 m s⁻¹ for the ice monomers. ⁸ This low critical velocity is also one of the obstacles in the planet formation in circumstellar disks. In this follow-up paper, we adopt v_crit ≃ 80 m s⁻¹ as a constraint on the growth of the dust aggregates.

Another limitation of Paper I is that the size of dust aggregates a_d and collisional velocity Δv show runaway growth in the collisional compression phase around S_t ∼ 1. However, this rapid growth would not be realistic if a more natural treatment of the internal density of the dust is considered (Section 2.2.1, see also Section 3).

Moreover, there is a critical process that may promote blanet formation. In paper I, we did not take into account the radial drift of the dust particles as the first approximation. The radial velocity of the dust v_r,d relative to the gas (Weidenschilling 1977; Tsukamoto et al. 2017) is v_r,d ∼ S_t η v_K , and $\eta \sim {({c}_{s}/{v}_{K})}^{2}$, where c_s is the gas isothermal sound velocity and v_K is the Keplerian rotational velocity. In the circumnuclear disk around an SMBH, initially S_t η ∼10⁻⁴–10⁻³. Then the drift time of the dust particle t_drift ∼ r/v_r,d ∼ 5–50 ${({M}_{\mathrm{BH}}/{10}^{7}{M}_{\odot })}^{1/2}{(r/1\mathrm{pc})}^{1/2}$ Myr. This is not negligibly small for the lifetime of the active galactic nucleus (AGN), i.e., 10⁷ − 10⁸ yr. In this paper, we investigate the effects of the radial drift of the dust particles.

The remainder of this paper is organized as follows. In Section 2, we describe the models for the dust evolution and its application to the circumnuclear region. In Section 3, we show the results of the models with and without the radial drift of the dust particles. In Section 4, we discuss how the maximum collisional velocity and the formation timescale of blanets depend on the parameters α and M_BH. We also discuss the expected mass of the blanets and their radial distribution. Finally, we summarize the results in Section 5.

2.1. The Region of "Blanet Formation"

Here we briefly summarize the concept of dust evolution around SMBHs, as discussed in Paper I. Figure 1 shows a schematic of the AGN and the circumnuclear disk. An SMBH (with a mass of 10⁶–10¹⁰ M_⊙) is surrounded by an accretion disk, which radiates enormous energy (the bolometric luminosity is ∼10⁴²–10⁴⁵ erg s⁻¹), mostly as ultraviolet and X-rays. The dust particles in the central r < r_sub are sublimated by the radiation from the accretion disk around the SMBH. The radius depends on the AGN luminosity:

$$\begin{eqnarray}&&{r}_{{\rm{sub}}}\simeq 1.3\,{\rm{pc}}{\left(\displaystyle \frac{{L}_{{\rm{UV}}}}{{10}^{46}{\rm{erg}}{{\rm{s}}}^{-1}}\right)}^{0.5}{\left(\displaystyle \frac{{T}_{{\rm{sub}}}}{1500{\rm{K}}}\right)}^{-2.8}{\left(\displaystyle \frac{{a}_{d}}{0.05\mu {\rm{m}}}\right)}^{-1/2},\end{eqnarray} \tag{ 1 }$$

where L_UV is the ultraviolet luminosity of the AGN, and a_d is the dust size (Barvainis 1987). The radiation forms conical ionized gas (narrow emission-line region) and also contributes to producing the outflows of the dusty gas and torus (Wada 2012; Izumi et al. 2018; Wada et al. 2018). In the mid-plane of the torus, cold, dense gas forms a thin disk, where icy dust particles can be present beyond the snowline r_snow (see Section 2.4).

Zoom In Zoom Out Reset image size

We introduce a free parameter $\alpha \equiv {v}_{t}^{2}/{c}_{s}^{2}$, where c_s is the gas sound velocity and v_t is the turbulent velocity, to represent strength of the kinematic viscosity due to the turbulence (e.g., Ormel & Cuzzi 2007, see also Equations (7) and (8)). In the circumnuclear disk in AGNs, the value of α is highly uncertain. ⁹ Therefore, we here treat α as a free parameter to check how it alters the results, especially the maximum collisional velocity between the dust aggregates and the onset of the gravitational instability of the dust disk.

In contrast to the dust coagulation process in the circumstellar disks (Weidenschilling 1977), the drag between dust particles and gas obeys the Epstein law. The aggregate's size (a_d ) is always much smaller than the mean free path of the gas, ${\lambda }_{g}\sim {10}^{12}\ \mathrm{cm}{({\sigma }_{\mathrm{mol}}/{10}^{-15}{\mathrm{cm}}^{2})}^{-1}{({n}_{\mathrm{mol}}/{10}^{3}{\mathrm{cm}}^{-3})}^{-1}$, where σ_mol and n_mol are the collisional cross-section and number density of the gas, respectively.

2.2. Evolution of Dust Aggregates in Each Stage

The model for the growth of dust particles here is based on the elementary processes found around stars. The evolution of dust particles is divided into four stages as described below.

2.2.1. Hit-and-stick Stage

If the dust aggregates grow through ballistic cluster-cluster aggregation (BCCA), the internal structure of the aggregate should be porous (i.e., the internal density is much smaller than the monomer's density: ρ_int ≪ ρ₀), and its fractal dimension is D ≃ 1.9 (Mukai et al. 1992; Okuzumi et al. 2009). This is called the hit-and-stick stage (Okuzumi et al. 2012; Kataoka et al. 2013), and the internal density is given by

$$\begin{eqnarray}&&{\rho }_{\mathrm{int}}={\rho }_{0}{\left(\displaystyle \frac{{m}_{d}}{{m}_{0}}\right)}^{1-3/D},\end{eqnarray} \tag{ 2 }$$

where m_d is the mass of the aggregate, and m₀ is the monomer's mass. We assume that m₀ = 10⁻¹⁵ g and ρ₀ = 1 g cm⁻³.

The growth rate for m_d is

$$\begin{eqnarray}&&\displaystyle \frac{{{dm}}_{d}}{{dt}}=\displaystyle \frac{2\sqrt{2\pi }\,{{\rm{\Sigma }}}_{d}\,{a}_{d}^{2}\,{\rm{\Delta }}v}{{H}_{d}},\end{eqnarray} \tag{ 3 }$$

where Σ_d is the surface density of the dust disk, Δv is collision velocity between the aggregates, and H_d is the scale height of the dust disk as given in (Youdin & Lithwick 2007; Tsukamoto et al. 2017)

$$\begin{eqnarray}&&{H}_{d}={\left(1+\displaystyle \frac{{S}_{t}}{\alpha }\displaystyle \frac{1+2{S}_{t}}{1+{S}_{t}}\right)}^{-1/2}{H}_{g},\end{eqnarray} \tag{ 4 }$$

where H_g = c_s /Ω_K is the scale height of the gas disk, and S_t is the Stoke parameter, i.e., the normalized stopping time is defined as S_t ≡ t_stop/t_L with the eddy turn over time t_L . Here we suppose that ${t}_{L}={{\rm{\Omega }}}_{K}^{-1}$ and

$$\begin{eqnarray}&&{S}_{t}=\displaystyle \frac{\pi {\rho }_{{\rm{int}}}\,{a}_{d}}{2{{\rm{\Sigma }}}_{g}},\end{eqnarray} \tag{ 5 }$$

where Σ_g is the gas surface density.

The collision velocity between aggregates Δv for S_t < 1 can be divided into two regimes (Ormel & Cuzzi 2007): regime I) ${t}_{s}\ll {t}_{\eta }={t}_{L}\,{{Re}}^{-1/2}$, and regime II) ${t}_{\eta }\ll {t}_{s}\ll {{\rm{\Omega }}}_{K}^{-1}$. Here the Reynolds number, ${R}_{e}\equiv \alpha {c}_{s}^{2}/({\nu }_{\mathrm{mol}}{{\rm{\Omega }}}_{K})$ with the molecular viscosity ${\nu }_{\mathrm{mol}}\sim \tfrac{1}{2}{c}_{s}{\lambda }_{g}$ is

$$\begin{eqnarray}&&{R}_{e}\approx 9\times {10}^{4}\alpha {\left(\frac{{M}_{{\rm{BH}}}}{{10}^{6}{M}_{\odot }}\right)}^{-1/2}{\left(\frac{r}{1{\rm{pc}}}\right)}^{3/2}\left(\frac{{c}_{s}}{1\,{\rm{km}}\,{{\rm{s}}}^{-1}}\right).\end{eqnarray} \tag{ 6 }$$

For the hit-and-stick stage, ${S}_{t}\ll {R}_{e}^{-1/2}$, then for the regime I,

$$\begin{eqnarray}&&{\rm{\Delta }}{v}_{I}\simeq \sqrt{\alpha }{c}_{s}{R}_{e}^{1/4}| {S}_{t,1}-{S}_{t,2}| ={C}_{I}\sqrt{\alpha }{c}_{s}{R}_{e}^{1/4}{S}_{t},\end{eqnarray} \tag{ 7 }$$

where S_t,1 and S_t,2 are Stokes numbers of the two colliding particles, and C_I is a constant of the order of unity (Sato et al. 2016). For regime II, on the other hand,

$$\begin{eqnarray}&&{\rm{\Delta }}{v}_{{II}}\simeq {v}_{L}\sqrt{{t}_{\mathrm{stop}}/{t}_{L}}\simeq \sqrt{\alpha {S}_{t}}\,{c}_{s}\end{eqnarray} \tag{ 8 }$$

where v_L is velocity of the largest eddy. We assume that Δv_I = Δv_II at the transition.

The size of dust aggregates determines how they interact with the gas. The dynamics of the aggregates is affected by their cross sections, which depend on their internal inhomogeneous structure. The radius of BCCA cluster a_BCCA consisted of N monomers (N = m_d /m₀) is given as a_BCCA ≃ N^0.5 a₀ for N ≫ 1 (Mukai et al. 1992; Wada et al. 2008, 2009), which was also confirmed through N-body simulations (Suyama et al. 2012). We then assume that

$$\begin{eqnarray}{a}_{d}\left({m}_{d}\right)=\left\{\begin{array}{lc}{\left(\displaystyle \frac{{m}_{d}}{{m}_{0}}\right)}^{1/2}{a}_{0} & \left(\mathrm{BCCA}\right)\\ {\left(\displaystyle \frac{3{m}_{d}}{4\pi {\rho }_{\mathrm{int}}}\right)}^{1/3} & \left(\mathrm{otherwise}\right).\end{array}\right.\end{eqnarray} \tag{ 9 }$$

2.2.2. Compression Stages

The hit-and-stick stage ends due to collisions between the aggregates (collisional compression), or due to their interaction with the ambient gas (gas-drag compression). In the collisional compression, the rolling energy E_roll, which is the energy required to rotate a particle around a connecting point by 90^°, is comparable to the impact energy, E_imp = m_d Δv²/4, between the two porous dust aggregates of the same mass, m_d . Beyond this point, the aggregates start to get compressed due to mutual collisions and interaction with the gas (i.e., the ram pressure).

According to Suyama et al. (2012), the internal density of the aggregated ρ_int,f formed by collisions between two equal-mass aggregates, with density ρ_int, is calculated for E_imp > 0.45 E_roll:

$$\begin{eqnarray}&&{\rho }_{{int},f}={\left({\rho }_{\mathrm{int}}^{4}+{\rho }_{f}^{4}\displaystyle \frac{{E}_{\mathrm{imp}}-0.45{E}_{\mathrm{roll}}}{0.15{{NE}}_{\mathrm{roll}}}\right)}^{1/4},\end{eqnarray} \tag{ 10 }$$

ρ_f is the fractal density of the dust aggregate: ${\rho }_{f}\approx {m}_{d}/(7.7{a}_{d}^{2.5})$, and E_roll = 4.37 × 10⁻⁹ erg.

Moreover, the fluffy dust aggregates can be compressed owing to the ram pressure of the ambient gas (Kataoka et al. 2013). The internal density of the aggregates that are compressed by the gas is given

$$\begin{eqnarray}&&{\rho }_{\mathrm{int},\mathrm{drag}}\simeq {\left(\displaystyle \frac{{a}_{0}^{3}}{{E}_{\mathrm{roll}}}{P}_{g}\right)}^{1/3}{\rho }_{0},\end{eqnarray} \tag{ 11 }$$

where the ram pressure for a dust aggregate is

$$\begin{eqnarray}&&{P}_{g}=\displaystyle \frac{{m}_{d}\,{v}_{d}}{\pi {a}_{d}^{2}\,{t}_{s}},\end{eqnarray} \tag{ 12 }$$

with the stopping time t_s = S_t /Ω_K (Kataoka et al. 2013).

As the aggregates become more massive (m_d > 10¹⁰ g), they start getting compressed owing to their self-gravity, and the internal density evolves as ${\rho }_{\mathrm{int}}\propto {({\rm{\Delta }}v)}^{3/5}\,{m}_{d}^{-1/5}$ (Okuzumi et al. 2012).

2.2.3.  N-body Stage

When S_t ≃ 1, kinematics of the aggregates is affected not only by the turbulence, but also by mutual interaction between the aggregates as an N-body system and by the gravitational interaction with the density fluctuation due to the turbulence. Then the collision velocity between the aggregates is determined by a balance between various heating and cooling processes as the N-body particles. According to Michikoshi & Kokubo (2016, 2017), we solve the following equation to get the equilibrium random velocity of the dust aggregates v_d ,

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{{{dv}}_{d}^{2}}{{dt}} & = & {\left(\displaystyle \frac{{{dv}}_{d}^{2}}{{dt}}\right)}_{\mathrm{grav}}+{\left(\displaystyle \frac{{{dv}}_{d}^{2}}{{dt}}\right)}_{\mathrm{turb},\mathrm{stir}}+{\left(\displaystyle \frac{{{dv}}_{d}^{2}}{{dt}}\right)}_{\mathrm{turb},\mathrm{grav}}\\ & & -{\left(\displaystyle \frac{{{dv}}_{d}^{2}}{{dt}}\right)}_{\mathrm{coll}}-{\left(\displaystyle \frac{{{dv}}_{d}^{2}}{{dt}}\right)}_{\mathrm{drag}}=0.\end{array}\end{eqnarray} \tag{ 13 }$$

The first three heating terms represent the gravitational scattering of the aggregates, stirring by the turbulence, and gravitational scattering by density fluctuation of the turbulence, respectively. The two cooling terms in Equation (13) represent the collisional damping and the gas drag. We assume the collision velocity Δv ≈ v_d at S_t = 1 and numerically solve Equation (13) in this stage.

2.2.4. Radial Drift of the Dust Particles

In Paper I, we ignored the radial drift of the dust particles in the disk. However, as mentioned in Section 1, this is not always obvious. Here, we solve the following governing equations for the dust evolution based on the assumption that the mass distribution of the dust particles at each orbit radius is singly peaked at a mass (Sato et al. 2016; Tsukamoto et al. 2017);

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{\partial }}{{\rm{\Sigma }}}_{d}}{{\rm{\partial }}t}+\displaystyle \frac{1}{r}\displaystyle \frac{{\rm{\partial }}}{{\rm{\partial }}r}({{rv}}_{r,d}\,{{\rm{\Sigma }}}_{d})=0,\end{eqnarray} \tag{ 14 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial {m}_{d}}{\partial t}+{v}_{r,d}\displaystyle \frac{\partial {m}_{d}}{\partial r}=\displaystyle \frac{{m}_{d}}{{t}_{\mathrm{coll}}}.\end{eqnarray} \tag{ 15 }$$

Here t_coll in Equation (15) is the collision time, and the source term is

$$\begin{eqnarray}&&\displaystyle \frac{{m}_{d}}{{t}_{{\rm{coll}}}}={m}_{d}(4\pi {a}_{d}^{2}\,{n}_{d}\,{\rm{\Delta }}v)=2\sqrt{2\pi }{a}_{d}^{2}\,{{\rm{\Sigma }}}_{d}\,{\rm{\Delta }}{{vH}}_{d}^{-1},\end{eqnarray} \tag{ 16 }$$

where n_d is the number density of the dust particles.

The dust particles have a radial velocity due to the drag with the ambient gas:

$$\begin{eqnarray}&&{v}_{r,d}=-\left(\displaystyle \frac{{v}_{r,g}}{1+{S}_{t}^{2}}+\displaystyle \frac{2{S}_{t}}{1+{S}_{t}^{2}}\,\eta \,{v}_{{\rm{K}}}\right),\end{eqnarray} \tag{ 17 }$$

where v_K is the Kepler velocity and η is a parameter that determines the sub-Keplerian motion of the gas, and the radial velocity of the gas v_r,g is given with the mass accretion rate $\dot{M}$:

$$\begin{eqnarray}&&{v}_{r,g}=-\displaystyle \frac{\dot{M}}{2\pi r{{\rm{\Sigma }}}_{g}},\end{eqnarray} \tag{ 18 }$$

where the mass accretion rate $\dot{M}$ is assumed to be using the Eddington mass accretion rate $\dot{M}={\gamma }_{\mathrm{Edd}}\,{\dot{M}}_{\mathrm{Edd}}$ with the Eddington ratio γ_Edd.

2.3. Gravitational Instability of the Dust Disk and Formation of Blanets

We investigate the gravitational instability (GI) of the disk consisting of dust aggregates with S_t > 1 using the Toomre's Q-value defined as

$$\begin{eqnarray}&&{Q}_{d}\equiv \displaystyle \frac{({v}_{d}/\sqrt{3}){{\rm{\Omega }}}_{K}}{3.36G{{\rm{\Sigma }}}_{d}}.\end{eqnarray} \tag{ 19 }$$

For the axisymmetric mode, Q_d < 1 is the necessary condition for GI, but the non-axisymmetric mode can develop for Q_d ≲ 2. In this case, spiral-like density enhancements are formed followed by fragmentation of the spirals (Michikoshi & Kokubo 2017), which leads formation of massive objects, i.e., blanets. The mass of blanets is estimated as

$$\begin{eqnarray}&&{M}_{{\rm{bl}}}\simeq {\lambda }_{{\rm{GI}}}^{2}{{\rm{\Sigma }}}_{d},\end{eqnarray} \tag{ 20 }$$

where the critical wavelength for GI is

$$\begin{eqnarray}&&{\lambda }_{{\rm{GI}}}=\displaystyle \frac{4{\pi }^{2}G\,{{\rm{\Sigma }}}_{d}}{{{\rm{\Omega }}}_{K}^{2}}.\end{eqnarray} \tag{ 21 }$$

2.4. Initial and Boundary Conditions

In all the models, the circumnuclear cold gas disk embedded in the geometrically thick torus (see Figure 1) is assumed to be gravitationally stable; the Toomre's Q-value, Q_g ≡ c_s Ω_K /π GΣ_g = 2. The gas sound velocity is assumed to be ${c}_{s}^{2}={k}_{{\rm{B}}}{T}_{g}/\mu {{\rm{m}}}_{H}$ with T_g = 100 K and μ = 2.0.

The Eddington ratio is assumed to be γ_Edd = 0.01. The AGN bolometric luminosity is then L_bol = 1.3 × 10⁴²erg s⁻¹ M_BH/(10⁶ M_⊙). The X-ray luminosity of the AGN, which is used to determine the snowline radius, L_X = 0.1L_bol. This can be attributed to the fact that the UV flux from the accretion disk is attenuated in the dense circumnuclear disk.

The snowline for a_d = 0.1 μm,

$$\begin{eqnarray}\begin{array}{rcl}{r}_{{\rm{snow}}} & \approx & 1.5\,{\rm{pc}}{\left(\displaystyle \frac{{L}_{X}}{1.3\times {10}^{41}{\rm{erg}}{{\rm{s}}}^{-1}}\right)}^{1/2}{\left(\displaystyle \frac{{T}_{{\rm{ice}}}}{170{\rm{K}}}\right)}^{-2.8}\\ & & \times {\left(\displaystyle \frac{{a}_{d}}{0.1\mu {\rm{m}}}\right)}^{-1/2}.\end{array}\end{eqnarray} \tag{ 22 }$$

Therefore, it is expected that the dust in the most part of the circumnuclear disk is icy. We assume T_ice = 170 K. The dust-to-gas mass ratio is assumed as f_dg = 0.01. We solve the governing equations (Section 2.2.4) between r = 0.1 pc and 200 pc with 600 grid cells.

Figure 2(a) shows a typical evolution of a dust aggregate at the snowline for M_BH = 10⁶ M_⊙ and α = 0.02. The internal density of the aggregate ρ_int is plotted as a function of its mass m_d . Initially, the internal density decreases monotonically from the monomer's initial density, i.e., ρ₀ = 1 g cm⁻³ to 4 × 10⁻⁶ g cm⁻³, as its mass increases from m_d ∼ 10⁻¹⁵ g to ∼10⁻⁵ g. At that instant, the size of the aggregate becomes ∼1 cm (see Figure 2(b)). After this hit-and-stick phase, the fluffy dust aggregates keep growing due to collisions in the turbulent gas motion until S_t ≃ 1. During this stage (m_d = 10⁻⁵ g to 10¹⁰ g), the aggregates are compressed mainly due to the gas drag (Section 2.2.2), and therefore ρ_int gradually increases. ¹⁰ After S_t ≳ 1, the aggregates behave as an N-body system under the effect of the turbulent fluctuation (Equation (13)). For m_d > 10¹⁰ g and S_t > 1, the aggregates are compressed due to their self-gravity. In the model shown in Figure 2, it becomes gravitationally unstable at t = 75 Myr (see also Section 4.1).

Zoom In Zoom Out Reset image size

Figure 2(b) plots the collisional velocity Δv of the aggregate and its size a_d as a function of m_d . The size a_d monotonically increases. In the compression stage (m_d > 10⁵ g), the increase of a_d slows down (see Equation (9)). Initially Δv is 20 cm s⁻¹, and it slightly decreases during the hit-and-stick stage. Then it turns to an increase phase until ∼57 m s⁻¹ around S_t = 1 during the compression stage. The size of the aggregate becomes a_d ∼ 10⁴ cm at the end of this stage. In this case, the aggregates are not compressed by their self-gravity (m_d < 10¹⁰ g) before the dust disk becomes GI. Note that the increase of Δv slows down at m_d ∼ 0.1 g, which corresponds to the transition between Δv_I and Δv_II (Equations (7) and (8)).

The growth time of an aggregate for ${\rm{\Delta }}v={\rm{\Delta }}{v}_{I}\,\simeq 1/2\sqrt{\alpha }{c}_{s}{R}_{e}^{1/4}{S}_{t}$ can be estimated as

$$\begin{eqnarray}&&\begin{array}{l}{t}_{{\rm{g}}{\rm{r}}{\rm{o}}{\rm{w}}}\,\equiv \,{\left(\displaystyle \frac{d{\rm{l}}{\rm{n}}{m}_{d}}{dt}\right)}^{-1}=\displaystyle \frac{4\sqrt{2\pi }}{3}\displaystyle \frac{{H}_{d}\,{\rho }_{{\rm{i}}{\rm{n}}{\rm{t}}}\,{a}_{d}}{{\rm{\Delta }}v\,{{\rm{\Sigma }}}_{d}}\\ \,\simeq \,\displaystyle \frac{16}{3}\sqrt{\displaystyle \frac{2}{\pi }}\displaystyle \frac{{H}_{g}}{\sqrt{\alpha }\,{R}_{e}^{1/4}\,{c}_{s}\,{f}_{dg}}\\ \,\simeq \,2.9\times {10}^{7}\,[{\rm{y}}{\rm{r}}]\,{\left(\frac{{c}_{s}}{1{\rm{k}}{\rm{m}}{{\rm{s}}}^{-1}}\right)}^{-1}{\left(\displaystyle \frac{{f}_{dg}}{0.01}\right)}^{-1}\left(\displaystyle \frac{{H}_{g}}{0.1\,{\rm{p}}{\rm{c}}}\right)\\ \,\times \,{\left(\displaystyle \frac{\alpha }{0.02}\right)}^{-1/2}{\left(\displaystyle \frac{{R}_{e}}{{10}^{4}}\right)}^{-1/4}.\end{array}\end{eqnarray} \tag{ 23 }$$

Figure 3 shows time evolution of ρ_int, a_d , S_t , and Δv for the same model shown in Figure 2. The hit-and-stick phase lasts for ∼10 Myr as expected by t_grow, and S_t becomes unity at t = 60 Myr. At this moment, the dust aggregate's size reaches ∼100 m (Figure 3(b)).

Zoom In Zoom Out Reset image size

Figure 3(b) shows that the growth of a_d is exponential, or slower in time, in contrast to the results in Paper I. This is a natural consequence of the evolution of the mass of the dust aggregates. The mass increase rate of the dust is

$$\begin{eqnarray}&&\displaystyle \frac{{{dm}}_{d}}{{dt}}\sim {n}_{d}\,{a}_{d}^{2}\,{\rm{\Delta }}v\sim \displaystyle \frac{{{\rm{\Sigma }}}_{d}}{{H}_{d}}\,{a}_{d}^{2}\,{\rm{\Delta }}v\propto {S}_{t}^{1/2}\,{a}_{d}^{2}\,{\rm{\Delta }}v.\end{eqnarray} \tag{ 24 }$$

Here we assume the dust layer is sedimented, i.e., its thickness H_d is scaled as ${H}_{d}\propto {S}_{t}^{-1/2}$. For ${\rm{\Delta }}v\propto {S}_{t}^{1/2}$ (Equation (8)), dm_d /dt ∝ m_d , therefore, m_d grows exponentially. The runaway growth seen in Paper I is caused by the assumption that the internal density of the dust aggregates stays porous (i.e., the fractal dimension is ∼2) through the evolution. In reality, when the compression by the ambient gas works, ρ_int is nearly constant (i.e., ${m}_{d}\propto {a}_{d}^{3}$), as shown in Figure 2, therefore S_t ∝ a_d . If the scale height of the dust disk is constant, then ${{dm}}_{d}/{dt}\propto {S}_{t}^{1/2}{a}_{d}^{2}\propto {m}_{d}^{5/6};$ therefore, the growth of the dust aggregate should be slower than $\exp (t)$.

Figure 3(c) shows that the collisional velocity Δv gradually decreases during the hit-and-stick stage, and it turns to rapid increase during the compression stage until S_t becomes unity at t = 56 Myr. In the N-body stage, the collisional velocity decreases from its maximum value, 57 m s⁻¹, and it becomes GI (i.e., Q_d ≤ 2) at t ≃ 75 Myr.

For comparison, the evolution of the model without the radial drift is shown in Figure 4. We found that the dust aggregates before S_t = 1 evolve almost the same way as that in the model with the radial drift (Figure 3). However, the time for GI is 136 Myr, in contrast to 71 Myr for the case with the advection.

Zoom In Zoom Out Reset image size

Figure 5 shows the radial distributions of Σ_d , H_d , Δv, the radial velocity of the dust and gas (v_r,d and v_r,g ), ρ_int, S_t , m_d , and a_d at 200 Myr in the same model shown in Figure 2. As Figure 5(a) shows, the dust is accumulated around r ∼ 2 − 3 pc, where v_r,d ≪ v_r,g (Figure 5(c)) and S_t turns from S_t < 1 to S_t > 1 (Figure 5(b)). From Figure 5(d), the dust aggregates evolve more rapidly in the inner region (r ≲ 3 pc), and the maximum size is ∼km. We call these objects as blanetesimals.

Zoom In Zoom Out Reset image size

4.1. Dependence on α and M_BH.

In the models with the radial drift of the dust aggregates, we investigated how the maximum velocity of the collision ${\rm{\Delta }}{v}_{\max }$ and the time for GI (τ_GI) depend on α and M_BH. In Figure 6, we plot ${\rm{\Delta }}{v}_{\max }$ and τ_GI as a function of α for M_BH = 10⁶ M_⊙ and 10⁷ M_⊙. It shows that ${\rm{\Delta }}{v}_{\max }$ depends on α, and not on M_BH. If ${\rm{\Delta }}{v}_{\max }\,\lesssim 80$ m s⁻¹ is necessary for collisional growth as numerical experiments suggested (Wada et al. 2009), then α should be ∼0.04 or smaller.

Zoom In Zoom Out Reset image size

The behavior of the dust growth (e.g., Figure 2) does not significantly depend on α and M_BH, but the timescale to reach S_t = 1 is different as shown in Figure 6(b). For example, for M_BH = 10⁶ M_⊙ and α = 0.02, it takes ∼60 Myr when S_t exceeds unity, whereas it is ∼225 Myr for M_BH = 10⁷ M_⊙ and α = 0.02. This implies that smaller BHs may preferentially host blanets within a lifetime AGNs (≲10⁸ yr). Figure 6(b) shows that, for M_BH = 10⁷ M_⊙, the blanetesimal disk may not become GI earlier than ∼150 yr for α < 0.04. For α > 0.05 or 0.06, GI does not occur at r = r_snow in the models with M_BH = 10⁶ M_⊙ or 10⁷ M_⊙, respectively.

4.2. Number and Mass of Blanets

In the final stage of the evolution, the blanetesimal disk can be gravitationally unstable, and it fragments into massive objects, i.e., blanets (see Section 2.3). Figure 7 shows the radial distribution of the mass and typical separation between blanets, λ_bl ≈ λ_GI (Equation (21)). Two models with the radial drift for M_BH = 10⁶ M_⊙ with α = 0.02 and M_BH = 10⁷ M_⊙ with α = 0.04 are shown. For comparison, a model without the radial drift is also shown (M_BH = 10⁶ M_⊙ and α = 0.02). The mass of blanets ranges from ≃ 20M_E at r = r_snow to ≃3000M_E at r ∼ 3.5 pc for M_BH = 10⁶ M_⊙, in contrast to the model without advection, which is M_bl ≃ 3M_E − 7M_E . For M_BH = 10⁷ M_⊙, M_bl ≳ 10⁴ M_E − 10⁵ M_E outside the snowline. However, this extraordinary massive blanet is unlikely, because it is comparable to the minimum mass of brown dwarfs (∼2 × 10⁴ M_E ). Therefore, the largest size of the blanets (R_bl) would be maximally ∼10 × Earth's radius at r ∼ 3 pc for M_BH = 10⁶ M_⊙, if the average internal density is similar to that of the Earth.

Zoom In Zoom Out Reset image size

Number of blanets is ∼ 8 ×10⁵ at r = r_snow for M_BH = 10⁶ M_⊙ with α = 0.02, provided that all the dust is converted to blanets through GI. From Figure 7, the average distance between blanets would be λ_bl ∼ 10⁻³ pc at r = r_snow. Therefore, the system of blanets does not resemble any known exoplanet systems, in a sense that the "planets" are isolated, dominant objects in their orbits.

Is the swarm of blanets the final form of the system? The collisional timescale t_coll for the blanets can be estimated as

$$\begin{eqnarray}&&{t}_{\mathrm{coll}}^{-1}=4\sqrt{\pi }\,{n}_{\mathrm{bl}}\,{\sigma }_{\mathrm{bl}}\left({r}_{\mathrm{coll}}^{2}+\displaystyle \frac{G\,{m}_{\mathrm{bl}}}{{\sigma }_{\mathrm{bl}}^{2}}{r}_{\mathrm{coll}}\right),\end{eqnarray} \tag{ 25 }$$

where where n_bl and σ_bl are the number density and velocity dispersion of the blanets, respectively. r_coll is the distance at the closest approach, which is r_coll ∼ 2R_bl (Binney & Tremaine 2008). The second term of Equation (25) represents the gravitational focusing, which enhances the collision rate. For ${n}_{\mathrm{bl}}\sim {\lambda }_{\mathrm{bl}}^{-3}\simeq {({10}^{-3}\mathrm{pc})}^{-3}$, M_bl ≃ 20M_E , R_bl ≃ 3R_E , and σ_bl ∼ Δv ≃ 30 m s⁻¹ (Figure 3(c)), the collisional time is t_coll ≈ 6 Gyr. This can be shorter, if three-body encounters between a close-binary and a blanet are considered. Therefore, the blanet system could dynamically evolve in the cosmological timescale, and it would be interesting study the evolution by direct N-body simulations.

4.3. Can Blanets Acquire a Massive Gas Envelope?

It would be interesting to investigate whether blanets can obtain an atmosphere. If the Hill radius r_H of a blanet is larger than the disk scale height, the blanets may form a gap, then the obtained mass of the atmosphere depends on r_H . The Hill radius of a blanet, ${r}_{H}={({M}_{\mathrm{bl}}/3{M}_{\mathrm{BH}})}^{1/3}\,r$ is

$$\begin{eqnarray}&&{r}_{H}\simeq {10}^{-4}\,\mathrm{pc}{\left(\displaystyle \frac{{M}_{\mathrm{bl}}}{1000{M}_{E}}\right)}^{1/3},\end{eqnarray} \tag{ 26 }$$

whereas the scale height of the gas disk at r = 1.5 pc for M_BH = 10⁶ M_⊙ is, ${H}_{g}\simeq {c}_{s}/{{\rm{\Omega }}}_{K}\simeq 2\times {10}^{-2}{M}_{\mathrm{BH},6}^{-1/2}$ pc, which is larger than r_H .

On the other hand, the Bondi radius is

$$\begin{eqnarray}&&{r}_{{\rm{B}}}\equiv \displaystyle \frac{{{GM}}_{\mathrm{bl}}}{{c}_{s}^{2}}\simeq 1.3\times {10}^{-5}\,\mathrm{pc}\left(\displaystyle \frac{{M}_{\mathrm{bl}}}{1000\,{M}_{E}}\right){\left(\displaystyle \frac{{c}_{s}}{1\mathrm{km}{{\rm{s}}}^{-1}}\right)}^{-2}.\end{eqnarray} \tag{ 27 }$$

Therefore, we suspect that the mass of the atmosphere is limited by the Bondi accretion, rather than by a gap formation. If this is the case, the envelope mass would be

$$\begin{eqnarray}&&{M}_{\mathrm{env}}\sim \displaystyle \frac{4\pi }{3}{r}_{{\rm{B}}}^{3}\,{\rho }_{g}\end{eqnarray} \tag{ 28 }$$

$$\begin{eqnarray}&&\approx {10}^{-7}{M}_{E}{\left(\displaystyle \frac{{M}_{\mathrm{bl}}}{1000{M}_{E}}\right)}^{3}{\left(\displaystyle \frac{{c}_{s}}{1\mathrm{km}{{\rm{s}}}^{-1}}\right)}^{-6}\left(\displaystyle \frac{{\rho }_{g}}{2\times {10}^{-21}\,{\rm{g}}\,{\mathrm{cm}}^{-3}}\right).\end{eqnarray} \tag{ 29 }$$

However, if the Bondi radius is filled with the ambient gas during the orbital motion of a blanet, the envelope mass could be maximally

$$\begin{eqnarray}&&{M}_{\mathrm{env}}\sim 2\pi \,{r}_{\mathrm{bl}}\cdot \,\pi {r}_{{\rm{B}}}^{2}\,{\rho }_{g}\sim 0.05{M}_{E},\end{eqnarray} \tag{ 30 }$$

for M_bl = 1000M_E at r_bl = 1.5 pc. The small M_env would suggest that the runway accretion of the gas to blanets and formation of massive "gas giants" are difficult, but it also depends on how quickly are the orbits of blanets filled with the gas. This is an interesting open question for future studies.

4.4. Can Other Mechanisms of Planet Formation Be Applicable?

In this paper, we focused on the evolution of the dust aggregates based on the coagulation theories of dust monomers and the gravitational instability. Other formation mechanisms of planetesimals have also been proposed and extensively discussed in the planet formation community. Among them, we here look over the pebble accretion, secular gravitational instability (secular GI), and the streaming instability as possible mechanisms of blanetesimal formation. More quantitative analysis in the circumnuclear environment would be interesting for future studies.

The pebble accretion is an accretion process of small solid bodies (i.e., pebbles) to massive seed objects (e.g., planets or planetesimals) under the effect of the gas drag and gravity (e.g., Ormel 2017; Lambrechts et al. 2019). In the present case, relatively massive aggregates could accumulate ambient smaller particles, then they could become more massive. However, the timescale conditions for the pebble accretion (Ormel & Liu 2018), i.e., t_settl < t_enc and t_stop < t_enc are not likely to be satisfied in the circumnuclear disk, where the settling time t_settl is the time needed for a particle to sediment to the massive objects, and the encounter time t_enc is the duration of the gravitational encounter time. In other words, the radial flux of pebbles due to the gas drag is too small in the region of S_t ≪ 1. Moreover, the relatively large turbulent motion of the gas (i.e., α ≳ 0.01) may prevent from the pebble accretion (Ormel & Liu 2018), in contrast to the circumstellar disk. Therefore, we do not expect that the pebble accretion is a major process as a formation mechanism of blanets.

The secular GI, which is the gravitational instability due to gas-dust friction, is another possible mechanism to form planets in the circumstellar disk (Youdin 2011). According to Takahashi & Inutsuka (2014), the condition for this instability is expressed as

$$\begin{eqnarray}&&{\rm{\Gamma }}\equiv \left(\displaystyle \frac{\alpha }{4\times {10}^{-5}}\right){\left(\displaystyle \frac{{f}_{{dg}}}{0.1}\right)}^{-2}\left(\displaystyle \frac{{Q}_{g}}{10}\right)\left(\displaystyle \frac{\eta }{0.01}\right)\lesssim 1.\end{eqnarray} \tag{ 31 }$$

In the present case, Γ ∼ 100 for α ≳ 0.01, f_dg = 0.01, η ∼ 10⁻⁴ and Q_g = 2, therefore, we do not expect the secular GI in the circumnuclear disk.

The streaming instability could be an effective mechanism to make condensations of dust particles, if the dust-to-gas mass ratio (f_dg ) is close to unity (e.g., Youdin & Goodman 2005). According to numerical simulations by Carrera et al. (2015), coagulation of particles driven by the streaming instability depends on f_dg and S_t . They found that the streaming instability occurs for f_dg ∼ 0.02 for S_t ∼ 0.1. In our case, the dust-to-gas mass ratio increases to f_dg ≃ 0.02 − 0.03 from the initial value (0.01) outside the snowline due to the radial drift of the dust aggregates. However, S_t ≫ 1 in the region (Figure 5(b)), therefore we do not expect that the streaming instability occurs in the circumnuclear disk.

In this follow-up paper of Wada et al. (2019) (Paper I), we theoretically investigated a process of dust evolution around an SMBH in the galactic center. We proposed that a new class of astronomical objects, blanets can be formed, provided that the standard scenario of planet formation is present in the circumnuclear disk.

Here, we investigated the physical conditions of the blanet formation outside the snowline (r_snow ∼ several parsecs) in more detail, especially considering the effect of the radial drift of the dust aggregates. We also improved the dust evolution model in Paper I in terms of the internal density evolution of the dust aggregates. We assumed the maximum collisional velocity for destruction, which was suggested by previous numerical simulations, as one of the necessary conditions for blanet formation. We found that a dimensionless parameter $\alpha ={v}_{t}^{2}/{c}_{s}^{2}$, where v_t is the turbulent velocity and c_s is the sound velocity, describing the effective angular momentum transfer due to the turbulent viscosity in the circumnuclear disk should be smaller than 0.04 for the black hole mass M_BH = 10⁶ M_⊙; otherwise, the dust aggregates could be destroyed due to collisions. The formation timescale of blanets τ_GI at r_snow is τ_GI ≃ 70–80 Myr for α = 0.01 − 0.04. The blanets (M_bl) are more massive for larger radii; they range from M_bl ∼ 20M_E − 3000M_E in r < 4 pc, in contrast to M_bl = 3 − 7M_E for the case without the radial drift.

The typical separations between the blanets, estimated from the wavelength of the gravitational instability, would be ∼10⁻³–10⁻² pc.

For M_BH ≥ 10⁷ M_⊙, the formation timescale is longer than ∼150 Myr for α ≤ 0.04. Although the GI of the blanetesimal disk takes place just outside the snowline (r = 4.7 pc), they should not be blanets because they are more massive than the typical brown dwarf mass (∼3 × 10⁴ M_E ). Note that AGNs are often heavily obscured with dense gas even for hard X-rays (Buchner et al. 2014) (N_H > 10²³ cm⁻², i.e., Compton-thick). If this is the case, the snowline is located at the inner region (e.g., r ∼ 2 − 3 pc), and as a result, blanets with M_bl ∼ 10M_E − 100M_E around M_BH = 10⁷ M_⊙ could be possible. Our results suggest that blanets could be formed around relatively low-luminosity AGNs during their lifetime (≲10⁸ yr). The system of blanets should be extraordinarily different from the standard Earth-type planets in the exoplanet systems. The blanets may acquire a rarefied atmosphere due to accretion of the gas in the circumnuclear disk (Section 4.3). The dynamical evolution of the swarm of blanets around an SMBH and whether they become more massive objects or are destroyed due to collisions may be an interesting subject for future studies (Section 4.2).

We would like to appreciate the anonymous referee's many valuable comments. The authors also thank Hidekazu Tanaka for many thoughtful comments. This work was supported by JSPS KAKENHI grant No. 18K18774.